Turbulent Output-Based Anisotropic Adaptation 


Michael A. Park* 

Jan- Renee Carlson* 

NASA Langley Research Center, Hampton, VA 23681 


Controlling discretization error is a remaining challenge for computational fluid dynam- 
ics simulation. Grid adaptation is applied to reduce estimated discretization error in drag 
or pressure integral output functions. To enable application to high 0(10' ) Reynolds num- 
ber turbulent flows, a hybrid approach is utilized that freezes the near-wall boundary layer 
grids and adapts the grid away from the no slip boundaries. The hybrid approach is not 
applicable to problems with under resolved initial boundary layer grids, but is a powerful 
technique for problems with important off-body anisotropic features. Supersonic nozzle 
plume, turbulent flat plate, and shock-boundary layer interaction examples are presented 
with comparisons to experimental measurements of pressure and velocity. Adapted grids 
are produced that resolve off-body features in locations that are not known a priori. 

I. Introduction 

Obtaining a suitable grid remains one of the most difficult parts of the entire Computational Fluid Dynam- 
ics (CFD) simulation process. Grids must resolve flow features to provide adequate control of discretization 
error and remain small enough to permit reasonable computation times. The AIAA Drag Prediction Work- 
shops 1 are well-documented examples of the difficulty of suitable grid generation even under the watchful 
eyes of experts. Mavriplis 2 shows that discretization error still dominates other error sources (e.g., turbu- 
lence model distance function calculation, thin-layer viscous approximation, level of artificial dissipation) 
even with extremely fine grids that are computable with today’s supercomputers. Chaffin and Pirzadelr 3 
detail the effort required to manually specify grid resolution for accurate three dimensional (3D) high-lift 
computations. The high degree of manual intervention required in these cases prohibits the use of auto- 
mated design tools. The manual interaction required to create a suitable grid is reduced with adaptive grid 
methods. Automated parameter studies 4 and design 5 are possible with an automated grid generation and 
adaptation process. Accurate final adapted solutions can be produced that are independent of coarse initial 
grids. 6-9 There are two main barriers to the use of grid adaptation in CFD simulation. The first is the 
ability of the adaptive indicator to produce an improved grid that will eventually yield a solution suitable 
for engineering application. The second is the robustness of grid adaptation mechanics. 

Adaptive grid methods are designed to specify a resolution request that equidistributes and minimizes 
error estimates. The improved resolution request is commonly based on local error estimates. 10-13 Uniformly 
reducing the errors associated with all local-error sources of the flow may not be optimal from an engineering 
context, where calculating an output functional may be of greater concern. An alternative method is to 
estimate the error in the calculation of a specified engineering output functional. 14-1 ' Output error indicators 
utilize the dual or adjoint solution of an output functional to account for the impact of local error as well as 
the transport of these local errors throughout the problem domain to improve the calculation of that output 
functional. This output-adaptive approach has been applied to finite-volume simulations in two dimensions 
(2D) 16 and 3D 18 with a parallel execution scheme. 19 

Adapting 3D grids in the neighborhood of curved boundaries to an anisotropic metric is a difficult 
problem. Most of the difficulty lies in simultaneously maintaining the boundary nodes on the geometry 
description and element validity. The projection of nodes to curved geometry has been applied as a post- 
processing step after surface refinement to a linear surface approximation. Directly moving the surface nodes 
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to the projected location without modification of the volume grid may result in invalid, inverted elements. 
The displacements required to project surface nodes can be distributed throughout the grid to prevent 
the invalidation of elements near surfaces. The volume displacement has been computed via Laplacian 
smoothing, linear springs, torsional springs, trusses, and linear elasticity. While these methods reduce the 
incidence of invalid element creation, they are not guaranteed to prevent invalid elements while projecting 
boundary nodes. Iterative methods have been developed to project the surface node to the “first problem 
plane” 20,21 (past the first element that becomes inverted during projection) with various combinations of 
edge swapping and node movement. However, grid regeneration 22 is reported to be necessary when this 
iterative method fails. 

This study avoids the issue of adapting high aspect ratio elements on curved boundaries by utilizing 
a hybrid approach which freezes the near- wall grid and adapts elements outside this zone. This simplifies 
the adaptation mechanics and increases adaptation robustness. This hybrid approach is not applicable to 
problems with under resolved initial boundary layer grids, because that portion of the grid is held constant 
during the simulation. However, the hybrid approach is a powerful technique for problems with important 
off-body anisotropic features, i.e., plumes, sonic boom propagation, and bow shocks. This study is an 
extension of 3D output-based adaptation to subsonic and supersonic turbulent flows. 


II. Flow and Adjoint Solver 


Fully Unstructured Navier-Stokes Three-Dimensional (FUN3D) is a suite of codes for finite-volume 
CFD. 23,24 The FUN3D website, http://fun3d.larc.nasa.gov, contains the user manual and an exten- 
sive list of references. FUN3D is able to solve incompressible, Euler, and Reynolds-averaged Navier-Stokes 
(R.ANS) flow equations, either tightly or loosely coupled to a turbulence model. Domain decomposition 
is employed to fully exploit the distributed architectures predominate in high-performance computing to 
increase problem size and reduce the execution time of the simulation process. The compressible flow solver 
employs an unstructured finite-volume method for conserved variables stored at each node in the grid, 

Q n = [p pu pv pw E i>} T , (1) 

where p is density, u, v , and w are velocity, and E is total energy per unit volume. The Spalart-Allmaras 25,26 
(S-A) turbulence model quantity v is also included. FUN3D can utilize general mixed element grids, but 
only tetrahedral grids are used in this study. All node-based variables in the grid Q are computed by driving 
the flow equation residual R to steady-state with an implicit point-iterative method. The present study 
employs the RANS equations coupled to the S-A turbulence model. 

After the steady-state flow solution is obtained, the discrete adjoint equations 27, 28 are solved to complete 
the dual problem. To derive these equations, the residual R is augmented with a source term S and the 
sensitivity of the output function with respect to this source term is sought. A potential application of this 
sensitivity is to model truncation error as a source term and determine the effect of this error on the output 
function /. A Lagrangian is formed from / with R included as a constraint, 

£(Q,S,X) = f(Q,S) + \ T R(Q,S). ( 2 ) 

differentiated with respect 

( 3 ) 


( 4 ) 


Each adjoint variable A £ lZ bN corresponds to a flow equation. The Lagrangian is 
to S, 

dC df fdQ\ T f df , T dR\ , T dR 
dS~dS + {ds) [dQ + dQ) + OS' 

The dependency on ^ can be eliminated by solving the adjoint equation, 



With the adjoint equation satisfied, the sensitivity of the output function can be expressed as, 

dC _ df T dR 
dS~ dS + X dS' 


( 5 ) 


The linear adjoint equations in Eq. (4) are solved with a dual-consistent time-marching method. 29,30 
The iterative method is constructed such that an adjoint solution exactly dual to a direct differentiation 
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approach is recovered at each time step, yielding a convergence rate which is asymptotically equivalent to 
that of the flow equations. 


III. Output-Based Adaptation 


Venditti 16 describes an output-based error estimation and adaptation scheme. This has been implemented 
in 3D 18 with a parallel 19 implementation. The adaptation intensity / is formed on an embedded grid, 


It = 


(A" - A L Y R{Q H )\ + \Rx{\ H ){Q H - Q L ) 


^L\T I 


( 6 ) 


where X H and X L are the adjoint solution and Q H and Q L are the flow solution interpolated to the embedded 
grid. The ( ) H interpolation is based on a quadratic fit of the original grid solution and its reconstructed 
gradients. The () L interpolation is linear. The flow Q and adjoint R\ residual operators are also on the 
embedded grid. The embedded grid is restricted to the original grid I n . The embedded grid required 
for the Venditti error estimator can be prohibitively large (eight times original grid) for large 3D problems, 
which effectively limits the problem size of the output adaptive method. A single grid error estimation 
procedure is under development. It has been applied to Euler problems, 31 but not evaluated for turbulent 
simulations. 

Venditti 16 provides a procedure to calculate a new grid spacing request from the adaptive indicator I. n 
and an error tolerance tolo. The adaptation indicator is summed to find the global indicator Iq = )) ~^I n ■ 
The ratio of the remaining error to a user-specified error tolerance toln is 


_ In 
tolo ' 


( 7 ) 


The ratio of the control volume indicator to an equal share of toln is 


N 

toln 


(8) 


where N is the number of control volumes. When a cost function does not have an intuitive error 
i.e., sonic boom surface pressure integrals, the toln is set to half In at each adaptive iteration. The 
isotropic element length h is computed with an estimate of the spacing on the original mesh h° and 
and local error ratios, 



tolerance, 
requested 
the global 

(9) 


where an exponent of lu = 0.20 is based on an a priori estimate of the spatial error convergence. The 
anisotropy of mesh elements is based on the the Mach Hessian, where the element size in the smallest 
spacing direction is dictated by the adjoint adaptation parameter. 

A metric tensor M is employed to define the desired multidimensional grid resolution because it is a 
natural way to express local interpolation error estimates of linear functions. 10,32-34 The symmetric positive 
definite matrix M in 3D has the diagonal decomposition 


M = X 



X T = xax t . 


(10) 


The eigenvectors X define an orthonormal basis with length specifications hi in this basis, 



This metric can be interpreted as an ellipsoid with major and minor axes of directions V, and lengths ft, . 10 
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To form the metric M, the Mach Hessian H is decomposed into eigenvalues A and eigenvectors X. The 
matrix M has the same eigenvectors as matrix H to keep the same principle directions and absolute value 
of the corresponding eigenvalues to make it symmetric positive definite. The adjoint adaptation parameter 
is incorporated into the Hessian framework by scaling the three eigenvalues so that the largest eigenvalue 
corresponds to the adjoint adaptation spacing requirement. The eigenvalues and corresponding eigenvectors 
are sorted so that |Ai| > |A 2 | > |A 3 |. Thus, after scaling, the largest eigenvalue is 

Ms) 2 ’ < i2 » 

where the specified element length h is computed by the output error estimation process, Eq. (9). 

Parallel grid mechanics described by Park and Darmofal 31-35 are used to modify the tetrahedral back- 
ground grid. The local elemental operators of node insertion, node movement, element swap, and element 
collapse are utilized to iteratively drive the edge lengths to unity in M . 

IV. Jet Plume Simulation 

The supersonic wind tunnel geometry is shown in Fig. 1. This case is chosen because it has a complex 
shock and expansion structures in locations that are not known a priori. These shocks and expansions also 
interact with viscous boundary and sheer layers in a way that makes manual specification of grid element 
size and stretching difficult even if the locations of these features are known. The nozzle 6 from Putnam 
and Capone 36 has a exit-to-throat area ratio of 1.722 and a boattail angle of 5 deg. The ratio of jet total 
pressure to free stream static pressure is 8.12. The computational domain is a quarter of the axisymmetric 
model and the blade sting is ignored. The free stream Mach number is 2.2. The Reynolds number (based on 
maximum nozzle diameter D = 15.24cm) is 1.86 x 10 6 . The S-A turbulence model is utilized. This particular 
turbulence model is not ideal for modeling the shear flows of the plume, but has the advantage of a verified 
adjoint implementation in FUN3D. The computation domain for the nozzle plenum begins at the screens of 
the wind tunnel model where a total pressure boundary condition is applied. 3 ' The outflow boundary state 
is extrapolated from the interior of the domain. Free stream boundary conditions are applied to the far-held 
boundaries. The symmetry planes of the configuration are treated with y- and ^-symmetry. 



Figure 1- Sketch of jet-engine exhaust-nozzle simulator. All dimensions in centimeters unless otherwise noted. 


Figure 1. Plume test geometry. 36 


Accurate predictions of near-held pressure signatures are obtained with the output-based adaptation 
algorithm, which reduces estimated error in an off-body pressure integral. The integral of quadratic pressure 
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deviation over a surface s in the domain is 


2 


/ 


i 


p- Poo 
Poo 


ds, 


(13) 


where A s is the area of the integration surface. This surface is a cylinder with radius of D where the pressure 
field was measured in the experiment. This focuses the adaptation on improving the calculation of pressure 
deviation from free stream near this surface. The original grid was generated with VGrid. 38 The problem 
was executed on 128 processors. The symmetry plane of the 3D grid is shown in Fig. 2. The original grid 
has 521,198 nodes and the adapted grid has 4,754,224 nodes. FUN3D utilizes a node-based algorithm, so 
the number of nodes is the number of control volumes. The boundary layer grid is frozen within 0.05cm 
(0.00328 nozzle diameters) of the body, which is well inside of the boundary layer and advancing layer grid 
topology. 



(a) Original grid. 



(b) Final adapted grid. 


Figure 2. Symmetry plane grid. 


Pressure on the original and adapted symmetry plane grids is shown in Fig. 3. The nose shock and ex- 
pansion are better defined and propagated a longer distance from the body as a result of the grid adaptation. 
The shock train in the plume is also better defined. This is particularly evident on the aft shock where the 
original grid began to coarsen. Details of the plume region are shown in Fig. 4 and Fig. 5. The adapted grid is 
aligned with the boundary layer and shear layer. It is also aligned to the shocks and expansions. The regions 
where shocks or expansions interact with the boundary layer no longer have preferred alignment direction 
and result in isotropic refinement. The pressure field becomes more refined as a result of the adaptation, see 
Fig. 5. 



(a) Original grid. 


(b) Final adapted grid. 


Figure 3. Pressure contours on the symmetry plane. 


The pressure signature is extracted along a line that is one nozzle diameter above the configuration center 
line, as shown in Fig. 6. The pressure signature on the final adapted grid is the solid line. The wind tunnel 
data is from Putnam and Capone 36 and contains measurement by two different probes. Previous compu- 
tational results by Castner 39 are shown as the dotted line. There is significant scatter in the experimental 
measurements, which may be due to irregularities on the tunnel walls or the blade sting. The pressure is 
not freestream ahead of the bow shock, which may indicate an offset in the experimental reference pressure. 
Castner did not model the fore body so the dotted line is constant forward of the boattail x/D = 7.6. The 
signature due to the boattail and jet plume is shown in Fig. 7. The scatter in the wind tunnel measurements 
is greatest forward of the boattail expansion. The current results (solid line) accurately predict the location 
of the boattail expansion and primary plume shock. The predicted pressure is slightly greater than the 
measurements in the boattail expansion. 
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(a) Original grid. (b) Final adapted grid. 


Figure 4. Symmetry plane grid nozzle exit detail. 



(a) Original grid. (b) Final adapted grid. 


Figure 5. Nozzle exit detail of pressure contours on the symmetry plane. 



x/D 


Figure 6. Entire pressure signature at a radius of one nozzle diameter above the center line. 



x/D 


Figure 7. Aft pressure signature at a radius of one nozzle diameter above the center line. 
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Putnam and Capone 36 provides total pressure boundary layer profiles at a number of circumferential 
locations (j>, where <j> = 0 is opposite the blade sting. These experimental boundary layer profiles were 
measured with zero boattail angle at the start of boattail for the simulated configuration. The measurements 
include a small range of Reynolds number between 1.55 and 2.06 million based on maximum nozzle diameter. 
The Reynolds number based on this x/D = 6.833 model station is 12.69 million. A comparison of the total 
pressure boundary layer profiles for the original and adapted grids is provided in Fig. 8. The total pressure of 
the supersonic portions of the CFD profile has been reduce by using the normal shock relations, because it is 
assumed that the flow would experience a normal shock at the leading edge of the boundary layer probe. The 
wind tunnel measurements from the wind tunnel report are provided at multiple <f> = 0, 45, and 90 to give 
an indication of the experimental uncertainty. The blade sting is located at (f> = 180 so the disturbed profile 
data at (f> = 135 and 180 are omitted. The adapted grid is frozen below 0.00328 nozzle diameters, which 
is well inside the boundary layer. The outer portion of the boundary layer is under resolved in the initial 
grid, resulting in thicker boundary layer than the wind tunnel measurements. The adapted grid predicts a 
boundary layer with approximately the same thickness as the wind tunnel measurements. The adapted grid 
CFD predicts higher total pressure recovery for y/D < 0.05. The turbulent velocity profile for a subsonic 
flat plate with adapted grid is examined in the next section to further investigate this over prediction. 



Figure 8. Total pressure profile at x=104cm. 

An existing inviscid output-adaptive cut-cell technique 31,40 is also applied to the configuration to give 
an indication of the viscous effects on the signature. This comparison is included because Euler methods 
are the primary tool for studying near-field sonic boom signatures. The primary plume shock location is 
further aft for the inviscid simulation, Fig. 9. The inviscid secondary shock is further forward. The turbulent 
calculation better predicts the location of the primary plume shock location. 



x/d 

Figure 9. Aft pressure signature at a radius of one nozzle diameter above the center line. 
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V. Turbulent Flat Plate 


A turbulent flat plate is simulated with the hybrid adaptive scheme to verify the accuracy of the technique 
for boundary-layer flows. The selection of this case is motivated by the over prediction of total pressure (as 
compared to wind tunnel measurements) for the boundary layer profile in Fig. 8. Skin friction drag on 
the plate is used as the output to drive adaptation. The problem was executed on 64 processors. Similar 
conditions were investigated by Andren and Park. 41 A 3D channel with a width of 0.1 is simulated to 
exercise the 3D adaptation mechanics. Symmetry boundary conditions are applied to both side walls of the 
channel. The flat plate has a length of two and the symmetry boundary condition is applied to the floor 
of the channel for a length of one upstream of the flat plate. The Mach number is 0.1 and the Reynolds 
number is 12,693,200. The initial cell height is 1.92 x 10 -6 with a 1.1 geometric growth rate and the grid is 
frozen below 4.8 x 10 -4 , which is equivalent to y + = 212. These conditions are chosen to correspond to the 
state of the boundary layer and grid of the plume case at the location that the velocity profile was extracted. 
The initial grid has 88,389 nodes and the adapted grid has 349,152 nodes. The original grid (constructed by 
subdividing structured grid hexahedra in to tetrahedra) and adapted grid are shown in Fig. 10. 
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(a) Original grid. 



Figure 10. Symmetry plane grid for turbulent flat plate. 


The velocity profile is extracted at the center of the plate, a distance of one from the plate leading edge. 
This profile is compared with experimental curve fits in Fig. 11. The experimental curve fit is given as, 


u 


+ 


y+ : y+ < 10; 

oil ln(y+) + 5 : y+ > 10 


(14) 


The output adapted CFD results approach the experimental curve fits for small y + and large y + until the 
edge of the boundary layer at y + = 1000. The comparison of CFD to the experimental curve fits is extremely 
favorable for the output-adaptive technique for these flow conditions. 



Figure 11. Flat-plate turbulent velocity profile. 
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VI. Shock-Boundary Layer Interaction 


The shock-boundary layer experiment is described in Humble, Scarano, and van Oudheusden. 42 This case 
is selected due to its inclusion in the AIAA Supersonic Shock-Boundary Layer Interaction (SBLI) Workshop. 
The computational domain is shown in Fig. 12, which consists of an extruded 2D converging-diverging 
nozzle exiting into the test section containing a wedge set at 8 deg incidence to the tunnel center line. The 
Reynolds number is 5,683 per cm. The total to freestream pressure ratio is 11.56 and the total to freestream 
temperature ratio is 2.01. A total pressure and total temperature boundary condition 37 is applied to the 
inlet face upstream of the wind tunnel throat. These total conditions are intended to result in an isentropic 
expansion to the Mach 2.25 test condition. This Mach number is combined with a stagnation pressure of 50.5 
kPa and stagnation temperature of 295 K to dimensionalize the velocity for comparison with Particle Image 
Velocimetry (PIV) measurements. No slip boundary conditions are applied to the walls of the tunnel and 
the wedge shock generator. The outflow boundary condition is extrapolated from the interior of the domain. 
Drag on the floor of the tunnel test section is used as the output to drive adaptation. The problem was 
executed on 256 processors. The physics of shock-boundary layer interaction is inherently time-dependent. 
However, a steady-state simulation was employed in this study to obtain mean values. 



Figure 12. Shock-boundary layer test geometry. 


The original grid and adapted grid velocity fields along the center plane of the tunnel are shown in Fig. 13. 
The u - velocity is given in units of meters per second. The shock and expansion structures are better defined 
for the adapted grid. The wedge shock generator wake is propagated a longer distance down stream. The 
shock-boundary layer interaction zone on the lower wall of the adapted grid has a larger region of reversed 
flow. The shock reflected off the lower boundary layer has increased in strength and is more clearly defined. 



(a) Initial grid. 



(b) Output adapted grid. 


Figure 13. Velocity for the center plane (m/s). 


The intersection of the original and adapted tetrahedral volume grids with the center plane of the wind 
tunnel are shown in Fig. 14. The original grid was generated with VGrid 38 and has an initial cell height of 
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0.005 cm with a growth rate of 1.15 and secondary rate acceleration of 0.08. The boundary layer grid was 
frozen within 0.3 cm of the no slip boundary conditions. The initial grid has 766,384 nodes and the adapted 
grid has 1,301,836 nodes. The outer portions of the incoming boundary layers on the floor and roof of the 
test section have been refined. The bow shock of the wedge shock generator is tightly refined to propagate 
the shock to the shock-boundary layer interaction zone. The boundary layer above and below the wedge 
shock generator are better resolved. The shock reflected off the floor is also refined with the expansion off 
the the rear of the wedge. The adapted floor boundary layer grid is more isotropic at and down stream of 
the shock interaction zone due to increased stream wise Mach variation. 



(b) Output adapted grid. 

Figure 14. Detail of the intersection of the volume grid and center plane. 


The adapted grid velocity and experimental PIV data are shown in Fig. 15. The u-velocity is given in 
units of meters per second. The PIV measurements are not available in the upper left of the figure. The 
incident boundary layer is slightly thicker for the PIV measurement. The adapted CFD predicts a larger 
region of reversed flow. The CFD incident shock location is slightly forward of the PIV measurement. This 
may be due to the steady-state approximation. 

The adapted grid CFD and PIV data are extracted at various x locations near the shock boundary layer 
interaction zone, see Fig. 16. The x = 260 and 280 stations are upstream of the boundary layer thickening 
due to the shock interaction. The boundary layer thickens at x = 300 and 310 with CFD predicting the 
reversed flow region to extend further upstream of the shock. The largest difference between the CFD 
prediction and the PIV measurement is at and down stream of x = 320, where the CFD predicts higher 
velocity than the PIV measurement. 
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(a) Experimental time-averaged PIV ^-velocity field. 



(b) Final adapted Steady-state CFD u - velocity field. 


Figure 15. Detail of the shock-boundary layer interaction zone (m/s). 



Figure 16. Velocity profiles at various stations for shock-boundary layer interaction zone (m/s). 
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VII. Conclusions 


Grid adaptation is driven by estimated discretization error in drag and pressure integral output functions. 
To enable application to high 0( 10 7 ) Reynolds number turbulent flows a hybrid approach is utilized that 
freezes the near-wall portion of boundary layer grids and adapts the grid away from the no slip boundaries. 
Supersonic nozzle plume, turbulent flat plate, and shock-boundary layer interaction examples are presented 
with comparisons to experimental measurements of pressure and velocity. The conditions of the turbulent flat 
plate matched the grid and Reynolds number at the location of the boundary layer profile of the plume case 
and verified that the adaptive approach provides the expected boundary layer velocity profile for subsonic 
flow. Adapted grids are produced that resolve off-body features in locations that are not known a priori. 
These features include shocks, expansions, boundary layers, and shear layers. The grids are anisotropically 
adapted to each of these features that has a significant impact on the specified output. The regions where 
these features interact lose a preferred anisotropic direction so they are refined isotropically. The hybrid 
approach is not applicable to problems with under-resolved initial boundary layer grids, but is a powerful 
technique for problems with important off-body anisotropic features. 
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